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Thermodynamic properties of chiral spin liquids are investigated for a variant of the Kitaev model defined on 
a decorated honeycomb lattice. Using the quantum Monte Carlo simulation, we find that the model exhibits a 
finite-temperature phase transition associated with the time reversal symmetry breaking, in both topologically 
trivial and nontrivial regions. While changing the exchange constants, the phase transition changes from contin¬ 
uous to discontinuous one, apparently correlated with the change in the excitations from Abelian to non-Abelian 
anyons. We show this coincidence by computing the topological quantities: the Chern number and the thermal 
Hall conductivity. In addition, we find, as a diagnostic of the chiral spin liquids, successive crossovers with 
multi-stage entropy release above the critical temperature, which indicates that the hierarchical fractionalization 
of a quantum spin occurs differently between the two regions. 
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Understanding of quantum spin liquids (QSLs) in magnets, 
where strong quantum fluctuations suppress magnetic order¬ 
ing even at the lowest temperature (T), has been one of the 
most challenging subjects in strongly correlated electron sys¬ 
tems nil. Among many possible realizations of QSLs, 
the chiral spin liquid (CSL), in which the time reversal sym¬ 
metry is broken, has attracted considerable attention in not 
only condensed matter physics but also quantum information. 
This is because it may have the excitations obeying the non- 
Abelian anyon statistics, which are utilized as the operators 
in topological quantum computing 0 . To explore this exotic 
state, quantum spin systems on geometrically frustrated lat¬ 
tices have been intensively studied thus far |0] . For instance, 
the possibility of CSLs has been studied theoretically in the 
Heisenberg model on a kagome lattice . Experimentally, 
a possible CSL was discussed for a metallic pyrochlore com¬ 
pound Pr 2 lr 207 isl]. 

Besides the analyses of geometrically-frustrated quantum 
magnets, a class of the models that have the exact CSL ^ound 
states has opened a new avenue in the study of CSLs SMIii. 
One of them was originally suggested by A. Kitaev yl] and 
studied in detail by H. Yao and S. Kivelson J^. This model is 
a variant of the honeycomb Kitaev model, which is defined on 
a decorated honeycomb lattice, composed by extending each 
honeycomb lattice site to a triangle. The exact solution of this 
model shows that the ground state accommodates two differ¬ 
ent types of the CSLs accompanied with Abelian and non- 
Abelian anyons as the elementary excitations. Interestingly, 
the non-Abelian CSL has topologically nontrivial Majorana 
fermion bands with a chiral edge mode. 

Although the exact solutions for the ground states serve as 
good references in the exploration of CSLs, it is a crucial is¬ 
sue how the CSLs behave against thermal fluctuations at flnite 
T. As the discrete chiral symmetry is broken in the ground 
state, one expects a phase transition to the CSL at a flnite 
T, even in two dimensions. An interesting issue is how the 
non-Abelian anyonic nature survives at flnite T, since the ro¬ 
bustness is crucial for the application in topological quantum 


computing iQ. It is also important to clarify the precursors 
of CSLs when approaching from high-T paramagnet, as they 
may provide a clue for searching CSLs in experiments. How¬ 
ever, flnite-T properties of CSLs remain elusive, mainly be¬ 
cause of the lack of well-controlled theoretical methods for 
the quantum spin systems in question. For instance, the con¬ 
ventional world-line quantum Monte Carlo (QMC) method, 
which is one of the most standard techniques for flnite T, does 
not work efficiently due to the negative sign problem. 

In this Letter, we elucidate the effect of thermal fluctua¬ 
tions on CSLs by unbiased numerical simulations. We ad¬ 
dress this problem in a representative model that has the exact 
CSL ground state, a variant of the Kitaev model on a decorated 
honeycomb lattice mentioned above. We compute the flnite-T 
properties of this model by adopting the QMC method formu¬ 
lated in the Majorana fermion representation, which does not 
suffer from the negative sign problem lll2l4l4l] . We And that 
both topologically trivial and nontrivial CSLs exhibit a phase 
transition at a flnite T associated with the time reversal sym¬ 
metry breaking. By changing the exchange parameters, the 
phase transition changes from second order to first order. This 
appears to be related with the change of topological nature 
of the chiral phases, that is, from the Abelian to non-Abelian 
CSLs, which we corroborate by calculating the Chern num¬ 
ber and the thermal Hall conductivity. In addition, we clarify 
how the CSLs develop from the high-T paramagnet while de¬ 
creasing T; there are several crossovers characterized by the 
multi-stage entropy release, reflecting hierarchical fractional¬ 
ization of a quantum spin, which appears in a different way 
between the Abelian and non-Abelian regions. 

We consider a variant of the Kitaev model defined on a dec¬ 
orated honeycomb lattice, whose Hamiltonian is given by 
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where crj represents the 7 component of the Pauli matrix 
describing an S = 1/2 spin located at a site j. The dec¬ 
orated honeycomb lattice consists of triangles connected by 
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bonds; the interaction is defined for nearest-neighbor (NN) 
7 bonds in the triangles, while J'^ for NN 7 ' bonds 

connecting the triangles, as shown in the inset of 

Fig.ffl For simplicity, we assume J = Jx = Jy = Jz and 
J' = = Jy = J' as in Ref. JH], and introduce the parame¬ 

ter a so that J = cos a and J' = sin a. 

The model in Eq. o is exactly solvable for the ground 
state |@]. There are two kinds of local conserved quantities: 
one is defined for each triangle t, Wt = Ylj^t 
other for each dodecagon h, Wh = Ylj^h where 7 ^ de¬ 
notes the bond not included int or h among three NN bonds 
at the site j. Note that Wt changes its sign by the time re¬ 
versal operation as it consists of the three products of Pauli 
matrices. The exact ground state is a CSL in which the time 
reversal symmetry is broken by a uniform alignment of Wt 
as well as Wh. Remarkably, the ground state accommodates 
two topologically different CSLs depending on a, as pre¬ 
sented in the bottom of Fig. [T] The critical point is located 
Sit a = ac = 7r/3: the ground state is a topologically trivial 
CSL for a > ac, whereas it becomes topologically nontrivial 
for a < ac. The latter phase has the non-Abelian anyons in 
the excitation. 

In order to compute thermodynamic properties of CSLs in 
this model, we adopt a QMC technique developed by the au¬ 
thors and their collaborator recently imii. The method is 
based on the Majorana fermion representation of t he q uan- 
tum spins via the Jordan-Wigner transformation 1191. [iSl - tlTn . 
In terms of the Majorana fermions, the model in Eq. o is 
written as 


{jk)x {jk)y {jk)z 

+ij'^ cjCk - ij'y yy cjCk - iJi yy iir^cjCk, (2) 

U^)y 


where j < k. The operator r]r = icjCk defined on each 
z and z' bond is regarded as a Z 2 variable taking ±1, be¬ 
cause it commutes with the Hamiltonian and r]‘^ = 1 (r is the 
bond index). Thus, the Hamiltonian describes the free Ma¬ 
jorana fermions coupled with thermally-fluctuating Z 2 vari¬ 
ables. This representation enables the QMC simulation with¬ 
out the negative sign problem. We carried out 40,000 MC 
steps for measurement after 10,000 MC steps for thermaliza- 
tion. Moreover, we used the parallel tempering algorithm to 
avoid the slowing down at low T iiii: we prepared 16 repli¬ 
cas in each simulation. We calculated the 61/^-site clusters up 
to L = 10 with the twisted boundary condition 1I12I1 . 

Figure [T] shows the phase diagram obtained by the QMC 
simulation. We find that the model in Eq. o exhibits a phase 
transition at a finite T, as expected for the discrete chiral sym¬ 
metry breaking in the CSL phases. The critical temperature 
Tc is determined by the specific heat Cy and the chirality n, as 
described below. In addition to the transition, we find several 
crossovers as shown in the phase diagram; we will return to 
this point later. 

QMC data for sX a/ t: = 0.3 and 0.4 are shown in 



FIG. 1: (color online). Finite-T phase diagram of the Kitaev model 
on a decorated honeycomb lattice. The ground-state phase diagram 
is also presented in the bottom of the figure. The lattice structure 
is depicted in the inset. The circles represent the phase transition 
temperature Tc. The deduced location of the tricritical point is also 
shown; for larger (smaller) a, the transition is continuous (discon¬ 
tinuous). The triangles, squares, diamonds, and inverted triangles 
represent the crossover temperatures, 7 "* 7 "** 7 "*^ ^nd Th, respec¬ 
tively. The solid and dotted lines in the large a region represent Tc 
and T**, respectively, determined by the MC simulation for the ef¬ 
fective model for J'/J 1. The solid, dashed, and dashed-dotted 
lines in the small a region represent Tc, Tl, and respectively, 
obtained from the effective model for J'/J 1. See the text for 
details. 


Figs. [2a) and|2d), respectively. There is a sharp peak that 
grows with increasing the system size, indicating the phase 
transition. We also show the data for the mean squares of the 
chirality defined as ^ Wt in Figs. [2b) and [2^)- 
This quantity develops rapidly with decreasing T at the peak 
temperature of Cy , which clearly indicates that the phase tran¬ 
sition is associated with the time reversal symmetry breaking. 

Interestingly, Tc changes continuously while changing a as 
shown in Fig.[TJ despite the topological change in the ground 
state at a = Q^c- We, however, find that the nature of the phase 
transition changes in the vicinity of this point. To see this, 
we calculate the energy histogram at several T near Tc. Fig¬ 
ure [2g) shows the data sX a/ tt = 0.3. The histogram shows 
a double peak structure, which indicates that the phase transi¬ 
tion is of first order. On the other hand, we cannot find such 
behavior sX a/ir = 0.4. Instead, we show that the finite-size 
scaling collapse of {n?) works well for the T = 6, 8, and 10 
clusters as shown in Fig.[2h), suggesting that the phase tran¬ 
sition is of second order sXa/i: = 0.4. In the scaling collapse, 
the optimization is carried out by the Bayesian scaling analy¬ 
sis |[l9ll and we obtain the critical exponents as 1/z^ = 1.09(9) 
and T] = 0.18(3). These exponents are close to those for the 
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FIG. 2: (color online). T dependences of (a) the specific heat, (b) the 
mean square of the chirality (c) the mean of Z 2 variables Wh per 
dodecagon, Wh = ^'^h ajn = 0.3. The corresponding 

data at a = 0.4 are shown in (d), (e), and (f). The vertical dashed 
(dashed-dotted) line indicates Tc (T**). (g) Energy histograms at 
several T in the vicinity of Tc at a/n = 0.3. (h) Scaling collapse for 
(/^^) with fjv — 1.09 and 77 = 0.18 at q/tt = 0.4. 
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FIG. 3: (color online). T dependences of the absolute values of (a) 
the Chern number iy and (b) the thermal Hall conductivity at 
a/n = 0.3. The vertical dashed line indicates Tc and the dotted 
line in (b) represents ttT 112. Corresponding data at a j it = 0.4 are 
shown in (c) and (d). 


2D Ising universality class, u = 1 and 77 = 1/4 i2^ . The 
results suggest the existence of the tricritical point between 
q^/tt = 0.3 and 0.4, which is close to ac, as shown in Fig.[Tl 
To examine the relation between the tricritical point and the 
topological nature of CSL phases below Tc, we calculate the 
topological quantities at finite T by QMC. For this purpose, 
we first obtain the information of Majorana fermion bands. In 
the ground state, the bilinear Majorana fermion Hamiltonian, 
which is given by Eq. (O with all 77 ^ = 1, is diagonalized as 


= Efe ciHkCk = En Efe Snk{‘2.flkfnk “ 1), where Ck 
is a set of the Fourier transforms of Cj in a unit cell and 
is the Bloch Hamiltonian. Here, we introduce fnk and 
as fermion operators for the n-th band. The summation 
is taken only for the states with positive eigenvalues Snk- In 
the finite-T calculations, we extend these definitions straight¬ 
forwardly by considering the supercell of 6T^ cluster ob¬ 
tained in the QMC simulation. 

First, we compute the Chern number, which 
is nonzero in the topologically nontrivial CSF 
ground state for a < ac m m. We ex¬ 
tend the definition to finite T as z/(T, {r]r}) = 

^En,kMEuk)E 


T.™ {Unk I *^0 I Urnh ) {Umk | | Ttnfc ) 


where V = T|, /p is the Fermi distribution function, and 
the one-particle energy is given by Enk = 2|5nfc|; \unk) is 
the eigenvector for for a given configuration of {r7r}; 
vi = dHkjdki (I = a, 6), where ka and are the projections 
of k onto the reciprocal primitive vectors (see the inset of 
Fig. m, and we assume these vectors are orthogonal lEI]. 
In the present calculation, we take = 10 and S = 10“^; 
we confirmed the convergence with respect to Lk and 6. 
FigureOa) shows the QMC data for z/(T) = (|z^(T, {77r})|), 
which is computed for 100 samples during 40,000 MC 
steps. At a/TT = 0.3, z^(T) decreases rapidly near Tc while 
increasing T, and the change becomes sharper for larger 
T, reflecting the flrst-order nature of the transition. This 
suggests a discontinuous change of u{T) from 1 to 0 at Tc in 
the thermodynamic limit. In contrast, z/(T) is always zero at 
aji: = 0.4, as shown in Fig.[3tc). 

Next, we compute the thermal Hall conductivity, 
which reflects the topological nature of the excita¬ 
tions, since the heat is carried by the itinerant Ma¬ 
jorana fermions {cj}. The thermal Hall conductivity 
kT^{T) is evaluated in a similar way to the Chern num¬ 
ber: n^^{T) = {\K^^{T,{r]r})\), whore K^^{T,{r]r}) = 

T rp 1 \ Tm ('^'nk\Va\Umk) {Umk\vb\'^nk) 

y Z^n,fc l^rnj^n ^ With 

C2{Enk) = /;,{£)} mm. Figurelb) 

shows the QMC results. Similarly to the Chern number z^(T), 
kT^{T) sharply changes near Tc at a/7r = 0.3. At low T, 
shows T-linear behavior with the quantized coefficient 


K, 

as = 7rT/12. On the other hand, is always zero at 
ajiT = 0.4, as shown in Fig. [3d). 

Thus, both topological quantities y and behave differ¬ 
ently at a /it = 0.3 and 0.4: they become nonzero below Tc 
for the former, while always zero for the latter. The results 
suggest that the topological nature of CSF changes between 
the both sides of the tricritical point at ^ o^c separating the 
continuous and discontinuous phase transitions. 

In addition to the change of the nature of phase transition 
and ordered phase, we And that the paramagnetic phase above 
Tc also shows distinct behavior for a < ac and a > ac'. sev¬ 
eral crossovers appear in a different way, as shown in Fig. [T] 
First of all, the system exhibits a crossover at T = T* ^ 1 
in both regions, which is almost constant while changing a. 
In addition, for a < a^wo obtain two crossovers at T = T^ 
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FIG. 4: (color online). T dependences of (a) the specific heat and (b) 
the entropy per site at a/7r = 0.05 and 0.4 in the L = 6 cluster. 


and Tl (T^ > Tl ). On the other hand, for a > ac, the system 
exhibits a crossover at T = T** below T*. These crossovers 
are signatures of thermal fractionalization of a quantum spin 
into Majorana fermions, as explained below. 

First, we discuss the crossovers for a > ac. FigureUa) dis¬ 
plays T dependence of at a/7r = 0.4. The data show three 
peaks: a broad peak appears at T = T** ^ 2x10 in be¬ 
tween Tc ^ 2x10“^ andT* ^ 1. The high-T crossover at T* 
comes from the kinetic motion of itinerant Majorana fermions 
{cj} in Eq. ^ on the inter-triangle J' bonds, which corre¬ 
sponds to the development of spin-spin correlations on the 
dimers. The entropy of ^ In 2 is released in this crossover, as 
shown in Fig.lUb). On the other hand, the low-T crossover at 
T** and the phase transition at Tc are caused by localized Ma¬ 
jorana fermions: T** corresponds to the coherent alignment 
of local conserved quantities {Wh} for dodecagons as shown 
in Fig.[2f), while Tc corresponds to the other local conserved 
quantities {Wt} for triangles [or equivalently the chiral or¬ 
der parameter shown in Fig.I^e)]. Correspondingly, the en¬ 
tropy is released by | In 2 and ^ In 2, as shown in Fig. lU^b). 
This multi-state entropy release is explained by considering 
the limit of a/ir 0.5 (J' ^ J), where an effective Hamil¬ 
tonian is obtained by the perturbation in terms of J/J' llz^ . 
Investigating the effective model by a classical Monte Carlo 
simulation, we obtained the asymptotic behavior of T** and 
Tc, both of which well agree with the QMC results for the 
model in Eq. (O, as shown in Fig. miH. Interestingly, T** 
appears to merge into Tc in the vicinity of the tricritical point; 
for a smaller a, the coherent alignments of {Wh} and {Wt} 
take place simultaneously in the first-order transition at Tc, as 
shown in Figs. [2b) and|2c). 

Next, let us discuss the crossovers for a < o^c- As shown in 
the data at a/7r = 0.05 in Fig.[4l the system shows four-stage 


entropy release in this region: ^ ln2 at T = T* ^ J ~ 1, 
|ln2 at T* - 10-\ |ln2 at Tj* - 2 x 10“^ and ^ln2 
at Tc < 10“^ (Tc is unreachable because of the severe slow¬ 
ing down in QMC). Similar to the case with J' ^ J, the 
crossover at T* comes from itinerant Majorana fermions, cor¬ 
responding to the spin-spin correlations are developed on the 
intra-triangle bonds. On the other hand, the two crossovers 
at Tl and T^ and the phase transition at Tc are explained by 
an effective model in the limit of a ^ 0 (J ^ J'), as fol¬ 
lows. The effective Hamiltonian is described by the pseudo 
spin T and Wt, which represent the four-fold degeneracy on 
each triangle at J' = 0 0. It has a similar form to the 
Kitaev model on a honeycomb lattice in the magnetic field 
/leff, whose ground state is topologically nontrivial and excita¬ 
tions are non-Abelian any on 13] . This consideration, together 
with the numerical results for the Kitaev model Q, suggests 
that the system shows two crossovers at T£ 0.048Jeff and 

.15 Jeff, where Jeff = J'/3, as a consequence of ther¬ 
mal fractionalization of pseudo spin r into two kinds of Majo¬ 
rana fermions. Moreover, Tc is also suggested as Tc oc J'^/ J 
from the effective model. The asymptotic behaviors are shown 
in Fig.dlby the dashed, dashed-dotted, and solid lines, which 
well agree with the QMC data for the model in Eq. (O. In¬ 
terestingly, this multi-stage entropy release indicates that the 
fractionalization of a quantum spin occurs differently in this 
non-Abelian region: the original spin cr is first fractionalized 
into a pseudospin r and Wt, and furthermore, r is fractional¬ 
ized into two Majorana fermions. 

In summary, we have clarified the finite-T properties of 
both topologically nontrivial and trivial CSLs realized in the 
Kitaev model on the decorated honeycomb lattice, by per¬ 
forming the QMC simulation in the Majorana fermion rep¬ 
resentation. We revealed that both CSLs exhibits the phase 
transition with chiral symmetry breaking at finite T, whose 
critical temperatures Tc seemingly meet with each other at 
the tricritical point. We corroborate this by computing the 
topological quantities such as the Chern number and the ther¬ 
mal Hall conductivity. We also find that the system exhibits 
several crossovers above Tc, reflecting the hierarchical frac¬ 
tionalization of a quantum spin into Majorana fermions. The 
present results promote understanding of both the CSL phases 
with Abelian and non-Abelian excitations at flnite T, and fur¬ 
thermore, the diagnostic of them in the high-T paramagnetic 
phase, which will be useful for the experimental exploration 
of CSLs in quantum magnets. 

We thank M. Udagawa for fruitful discussion. This work 
is supported by Grant-in-Aid for Scientiflc Research, the 
Strategic Programs for Innovative Research (SPIRE), MEXT, 
and the Computational Materials Science Initiative (CMSI), 
Japan. Parts of the numerical calculations are performed in 
the supercomputing systems in ISSP, the University of Tokyo. 


[1] P. W. Anderson, Mater. Res. Bull. 8, 153 (1973). 














5 


[2] L. Balents, Nature 464, 199 (2010). 

[3] A. Kitaev, Ann. Phys. 321, 2 (2006). 

[4] V. Kalmeyer and R. B. Laughlin, Phys. Rev. Lett. 59, 2095 
(1987). 

[5] K. Yang, L. K. Warman, and S. M. Girvin, Phys. Rev. Lett. 70, 
2641 (1993). 

[ 6 ] S. Gong, W. Zhu, and D. N. Sheng, Sci. Rep. 4, 6317 (2014). 

[7] B. Bauer, L. Cincio, B. P. Keller, M. Dolfi, G. Vidal, S. Trebst, 
and A W. W. Ludwig, Nat. Commun. 5, 5137 (2014). 

[ 8 ] Y. Machida, S. Nakatsuji, S. Onoda, T. Tayama, and T. Sakak- 
ibara. Nature 463, 210 (2010). 

[9] H. Yao and S. A. Kivelson, Phys. Rev. Lett. 99, 247203 (2007). 

[10] D. F. Schroeter, E. Kapit, R. Thomale, and M. Greiter, Phys. 
Rev. Lett. 99, 097202 (2007). 

[11] S. B. Chung, H. Yao, T. L. Hughes, and E.-A. Kim, Phys. Rev. 
B 81, 060403 (2010). 

[12] J. Nasu, M. Udagawa, and Y. Motome, Phys. Rev. Lett. 113, 
197205 (2014). 

[13] J. Nasu, M. Udagawa, and Y. Motome, J. Phys.: Conf. Ser. 592, 
012115 (2015). 

[14] J. Nasu, M. Udagawa, and Y. Motome, arXiv: 1504.1259 

[15] H.-D. Chen and J. Hu, Phys. Rev. B 76, 193101 (2007). 

[16] X.-Y. Peng, G.-M. Zhang, and T. Xiang, Phys. Rev. Lett. 98, 


087204 (2007). 

[17] H.-D. Chen, and Z. Nussinov, J. Phys. A Math. Theor. 41, 
075001 (2008). 

[18] K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 
(1996). 

[19] K. Harada, Phys. Rev. E 84, 056704 (2011). 

[20] The slight deviation of 77 may be due to the limited system sizes 
used in the scaling. Indeed, we obtained better agreement for 
larger system sizes for the effective model in the isolated dimer 
limit, J' ^ J. The details will be published elsewhere. 

[21] X.-P. Shi, Y. Chen, and J. Q. You, Phys. Rev. B 82, 174412 

( 2010 ). 

[22] J. M. Luttinger, Phys. Rev. 135, A1505 (1964). 

[23] C. L. Kane and M. P. A. Fisher, Phys. Rev. B 55, 15832 (1997). 

[24] A. Cappelli, M. Huerta, and G. R. Zemba, Nucl. Phys. B 636, 
568 (2001). 

[25] T. Qin, Q. Niu, and J. Shi, Phys. Rev. Lett. 107, 236601 (2011). 

[26] H. Sumiyoshi and S. Fujimoto, J. Phys. Soc. Jpn. 82, 023602 
(2013). 

[27] S. Dusuel, K. P. Schmidt, J. Vidal, and R. L. Zaffino, Phys. Rev. 
B 78, 125102 (2008). 

[28] J. Nasu and Y. Motome, unpublished. 


